REAL(8) FUNCTION objective_u_fe_SPAIN(x,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::bpval,wval,tauval,ct_global
REAL(8)::p0,m0,h0,q,b,objective_u_fe
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval
REAL(8)::aux_lhs,balt,aux0alt,minc,maxc,maxc2,maxc3
REAL(8)::term1,taxbar,bound1,bound2,new_value,old_value
REAL(8)::old_valuemax,sign_new,sign_old,sign_two,ctemp2,difc

REAL(8)::objective_u_feold,gn_global1old,x_global1old,ct_global1old 
INTEGER(4)::i_num_star,i_num_star2,index_optmax,i_num
REAL(8)::q2,aux02,p02,gnval2 ,temp1,temp2
INTEGER,PARAMETER::num=100
REAL(8)::ctemp(num),tolval

INTEGER(4)::i,jj,im,iax,feas,error,niter2,ftype,try,index_opt
LOGICAL::flag00
EXTERNAL q_fun

objective_u_feold = -1d10
i_num_star=0
tolval = 0d0

b0 = b_initial2
a0 = DEXP(a_initial2)

b = x
q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0           !p*gn - T
aux0 = aux0 + endowment
aux0 = aux0*(1d0-defaultgrid(i_default_global2))

ctval= a0+aux0
ct_global = ctval   

IF (ctval.LE.0d0) THEN
    objective_u_fe_SPAIN= huge_n +1d3*(0d0-ctval)           
    if (choice_glob.EQ.1) THEN
        gn_global1 = 0d0
        x_global1  = 0d0
        ct_global1 = 0d0
    ENDIF
    RETURN
ENDIF

!case 1: frictionless: no binding tax cap
!new with lumpsum taxes
IF (taxrate.GT.1d0) THEN
    try   = 100
ELSE
    try = 1
ENDIF

IF (chi.EQ.0d0) THEN
    cnval=1d0
    gnval = 0d0
ELSE
    term1 = (1d0-chi)/chi*(1d0-omegac)
    term1 = term1**(1d0/sigma)
    cnval = term1/(1d0+term1)
    gnval = 1d0-cnval
ENDIF


ctval = a0+aux0
hval  = 1d0


p0   = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
wval = alpha*p0
tauval = p0*gnval-aux0

2046 taxbar = taxrate*(a0+in*p0)

2047 flag00 = ((p0.LE.0d0).OR.(wval.LT.wbar).OR.(gnval.GT.1d0).OR.(ctval.LE.0d0)) 

IF (.NOT. flag00) THEN
    pres  = p0
    hres  = hval
    ctres = ctval
    taures= tauval
    wwres = wval
    gnres = gnval
    feas  = 1d0
ELSE
    ctres   = 0d0    
    feas    = 0d0
    taures  = tauval
    hres    = 0d0
    pres    = p0
    gnres   = gnval
    wwres   = wval
ENDIF

!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(wwres.LT.wbar)*(wbar-wwres)+(ctres.LE.0d0)*(0d0-ctres))
!without lumpsum taxes
IF (taxrate.LE.1d0) penalty = penalty +1d6*(taures.GT.taxbar)*(taures-taxbar)
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.1d0)*(gnres-1d0)+(gnres.LT.0d0)*(0d0-gnres))

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres  = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres  = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

cres = cres - (cnres.LE.1d-6)*(-1d6)
cres = cres - (ctres.LE.1d-6)*(-1d6)

IF (sigma.EQ.1d0) THEN
    utilc = DLOG(MAX(cres,1d-8))
ELSE 
    utilc = MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma) 
ENDIF

IF (sigmag.EQ.1d0) THEN 
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n +penalty 

objective_u_fe = wres

if (choice_glob.EQ.1) THEN
    gn_global1 = gnres
    x_global1  = hres
    ct_global1 = ctres
ENDIF

IF (try.EQ.1) THEN
    
    IF (objective_u_fe.GT.objective_u_feold) THEN
        objective_u_feold= objective_u_fe
        if (choice_glob.EQ.1) THEN
            gn_global1old = gn_global1    
            x_global1old  = x_global1 
            ct_global1old = ct_global1 
        ENDIF
    ENDIF

    !case 2: binding tax rate
    !pn*gn > 0
    
    temp1=-aux0/taxrate-a0
    IF (temp1.GT.0d0) THEN
        maxc = ((((1d0-omegac)/omegac)/temp1)**(1d0/(1d0+mu)))*ctval        !b2
        minc = 1d0-in*taxrate                                               !b6
    ELSE
        temp2 = temp1*taxrate/(in*taxrate-1d0)
        maxc  = (((1d0-omegac)/omegac/temp2)**(1d0/(1d0+mu)))*ctval         !b3
        maxc  = MIN(maxc,1d0-in*taxrate)                                    !b5
        maxc3=maxc
        IF (wbar.EQ.0d0) THEN
            minc  = 1d-5                                                    !b4
        ELSE
            minc  = 1d0-in*taxrate-(aux0+taxrate*a0)*alpha/wbar             !b4
        ENDIF   

    ENDIF    
    IF (wbar.EQ.0d0) THEN
        maxc2  = 1d0             
    ELSE                
        !alpha*pn >= wbar
        maxc2 = ((1d0-omegac)/omegac*alpha/wbar)**(1d0/(1d0+mu))*ctval      !b1
    ENDIF

    maxc2 = MIN(maxc,maxc2)
    maxc  = maxc2

    ftype = 402
    IF (i_num_star.EQ.0) old_value = AFunction(ftype,a0,aux0,ct_global,minc)
    old_valuemax = 10d+33 

    index_opt = 0
    index_optmax = 0
    difc = (maxc-minc)/(REAL(num)-1d0)

    150 DO i_num=1+i_num_star,num
        ftype=402
        ctemp(i_num) = minc+(REAL(i_num)-1d0)*difc     
        new_value = AFunction(ftype,a0,aux0,ct_global,ctemp(i_num))

        sign_new = SIGN(1d0,new_value)
        sign_old = SIGN(1d0,old_value)
        sign_two = SIGN(1d0,sign_new*sign_old)

        if (DABS(new_value) <= DABS(old_valuemax)) then
            index_optmax = i_num
            old_valuemax = new_value        
        end if      
        
        i_num_star = i_num
        if (sign_two.EQ.-1d0) then
            index_opt = i_num
            old_value = new_value
            !EXIT
        GOTO 222
        ENDIF
        old_value = new_value
    ENDDO       

    222 IF (index_opt.EQ.0) THEN
        bound1=ctemp(MAX(index_optmax-1,1))
        bound2=ctemp(MIN(index_optmax+1,num))
    ELSE
        bound1=ctemp(MAX(index_opt-1,1))
        bound2=ctemp(MIN(index_opt,num))        
    ENDIF

    IF (i_num_star.EQ.num)   try = 2
    
    cnval=BrentRoots(ftype,a0,aux0,ct_global,bound1,bound2,1d-10,1000,fval,niter2,error)  
    gnval = 1d0-cnval
   
    IF (error.EQ.1) THEN 
        gnval = -10d0
        IF (i_num_star.LT.num)  go to 150
        GOTO 2046
    ELSEIF (error.EQ.2) THEN
        PRINT*,'Problem FE: More Iterations!'
        PAUSE
        STOP
    ENDIF

    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    wval   = alpha*p0
    tauval = taxrate*(a0+in*p0)
    
    GO TO 2046

ELSE
    IF (objective_u_feold.GT.objective_u_fe) THEN
        objective_u_fe= objective_u_feold
        if (choice_glob.EQ.1) THEN
            gn_global1 = gn_global1old    
            x_global1  = x_global1old 
            ct_global1 = ct_global1old 
        ENDIF
    ENDIF
ENDIF

objective_u_fe_SPAIN = objective_u_fe

end FUNCTION objective_u_fe_SPAIN
!***************************************************************************
REAL(8) FUNCTION objective_u_unemp_SPAIN(x,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::lambda1,pp1,bpval,wval,tauval,objective_u_unemp
REAL(8)::p0,m0,h0,q,b,x_global2old,objective_u_unempold 
INTEGER,PARAMETER::num=100
REAL(8)::ctemp(num),aahat,aa1    
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval,gn_global2old,ct_global,ct_global2old
REAL(8)::aux_lhs,balt,qalt,aux0alt,minc,maxc,temp1,taxbar,netbc,bound1,bound2,new_value,old_value,old_valuemax
INTEGER(4)::i,jj,iax,feas ,error,niter2,try,ftype,i_num,index_opt,i_num_star,i_num_star2,index_optmax,root
real(8)::sign_new,sign_old,sign_two,ctemp2,difc,b_total,b_bailout
LOGICAL::flag00
EXTERNAL q_fun

objective_u_unempold = -1d10
i_num_star2 = 0
root=0

b0 = b_initial2
a0 = DEXP(a_initial2)
b = x
b_total = b 
q  = q_fun(b_total,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0 + endowment
aux0 = aux0*(1d0-defaultgrid(i_default_global2))

aux_lhs   = aux0+taxrate*(a0+in*wbar/alpha)                          !always smaller than cT
!aux_lhs < 0 follows from  pn*gn = aux0+tax < 0 
!                          alpha*p*h**alpha=wbar*h < wbar

IF (choice_nowbar.EQ.1) aux_lhs = 100d0     !no aux_lhs condition for wbar=0

ctval  = a0+aux0
ct_global = ctval

IF (choice_fe.EQ.1) THEN
    hres = 0d0
    gnres= 0d0
    wres = huge_n
    flag00=(hres.LE.0d0)
    try  = 1
    GO TO 2033
ELSEIF (choice_nowbar.EQ.1) THEN
    objective_u_unemp_SPAIN= huge_n
    if (choice_glob.EQ.1) THEN
        gn_global2= 0d0    
        x_global2 = 0d0
        ct_global2= 0d0
    ENDIF
    RETURN
ELSEIF  ((aux_lhs.LE.0d0).AND.(taxrate.LT.1d0)) THEN
    balt = b+1d-6
    q    = q_fun(balt,a_initial2) 
    aux0alt = -q*(balt-(1d0-lambda)*b0) + payoff*b0            
    aux0alt = aux0alt*(1d0-defaultgrid(i_default_global2))
    objective_u_unemp_SPAIN= huge_n+1d3*(aux0alt-aux0)             
    if (choice_glob.EQ.1) THEN
        gn_global2 = 0d0    
        x_global2  = 0d0
        ct_global2 = 0d0
    ENDIF
    RETURN    
ENDIF

!case 1: h<1 and nonbinding tax cap
try = 1

ftype = 101
minc  = 1d-3
!Alternative: minc = (alpha/wbar*(1d0-omegac)/omegac)**(1d0/(1d0+mu))*ctval                       !b8
maxc  = (wbar/alpha*omegac/(1d0-omegac)/(ctval**(1d0+mu)))**(-alpha/(mu*(alpha-1d0)+1d0+mu))      !b7

old_value = AFunction(ftype,a0,aux0,ct_global,minc)
old_valuemax = 10d+33 

i_num_star=0
index_opt = 0
index_optmax = 0
difc = (maxc-minc)/(REAL(num)-1d0)

454 DO i_num=1+i_num_star,num
    ftype=101
    ctemp(i_num) = minc+(REAL(i_num)-1d0)*difc     
    new_value = AFunction(ftype,a0,aux0,ct_global,ctemp(i_num))

    sign_new = SIGN(1d0,new_value)
    sign_old = SIGN(1d0,old_value)
    sign_two = SIGN(1d0,sign_new*sign_old)

    if (DABS(new_value) <= DABS(old_valuemax)) then
        index_optmax = i_num
        old_valuemax = new_value        
    end if      
    
    i_num_star = i_num
    if (sign_two.EQ.-1d0) then
        index_opt = i_num
        old_value = new_value
        !EXIT
    GOTO 555
    ENDIF
    old_value = new_value
ENDDO       

555 IF (index_opt.EQ.0) THEN
    bound1=ctemp(MAX(index_optmax-1,1))
    bound2=ctemp(MIN(index_optmax+1,num))
ELSE
    bound1=ctemp(MAX(index_opt-1,1))
    bound2=ctemp(MIN(index_opt,num))        
ENDIF

cnval=BrentRoots(ftype,a0,aux0,ct_global,bound1,bound2,1d-10,1000,fval,niter2,error)  


IF ((error.EQ.1)) THEN 
    cnval = -10d0
    IF (i_num_star.LT.num)  go to 454
ELSEIF (error.EQ.2) THEN
    PRINT*,'Problem: More Iterations!'
    PAUSE
    STOP
ENDIF

p0 = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
hval = (wbar/alpha/p0)**(1d0/(alpha-1d0))

1013 gnval  = hval**alpha-cnval
tauval = p0*gnval-aux0   

1015 taxbar = taxrate*(a0+in*p0*hval**alpha)
    
IF (chi.EQ.0d0) THEN
    flag00 = ((gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0)) 
ELSE
    flag00 = ((gnval.LT.0d0).OR.(gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0)) 
ENDIF

IF (.NOT.flag00) THEN
    pres  = p0
    hres  = hval
    gnres = gnval
    ctres = ctval
    taures= tauval
    wwres = wbar
    feas  = 1d0
ELSE
    ctres = 0d0 
    feas  = 0d0
    hres  = 0d0
    gnres = 0d0
    taures= tauval
    pres  = p0
    wwres = wbar    
ENDIF
  
!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(ctres.LE.0d0)*(0d0-ctres))
!without lumpsum taxes
IF (taxrate.LE.1d0) penalty = penalty+1d3*((taures.GT.taxbar)*(taures-taxbar))
penalty = penalty+1d9*((b_total.LT.bmin)*(bmin-b_total))+1d2*((b_total.LT.bmin))
penalty = penalty +1d3*((gnres.GE.hres**alpha)*(gnres-hres**alpha)+(gnres.LT.0d0)*(0d0-gnres))
penalty = penalty +1d3*((hres.GT.1d0)*(hres-1d0)+(hres.LE.0d0)*(0d0-hres))

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

IF ((cnres.LE.1d-6).OR.(ctres.LE.1d-6)) cres = 0d0

IF (alphac.EQ.1d0) THEN
    aahat  = 1d0
    aa1    = 1d0
ELSE
    aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
    aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))
ENDIF

    
IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n + penalty 

objective_u_unemp = wres
if (choice_glob.EQ.1) THEN
    gn_global2 = gnres
    x_global2 = hres
    ct_global2 = ctres
ENDIF

        
IF (try.EQ.1) THEN
    IF (root.EQ.0) THEN
        x_global2old  = x_global2
        gn_global2old = gn_global2
        ct_global2old = ct_global2
    
        objective_u_unempold = objective_u_unemp
        root=1
    ELSEIF (root.GT.0) THEN
        IF (objective_u_unemp.GT.objective_u_unempold) THEN
            objective_u_unempold= objective_u_unemp
        
            if (choice_glob.EQ.1) THEN
                gn_global2old = gn_global2    
                x_global2old  = x_global2 
                ct_global2old = ct_global2
            ENDIF
        ENDIF
        root=root+1
    ENDIF
    
    IF (i_num_star.LT.num)     GO TO 454
ENDIF

2033 IF (try.EQ.1) THEN
    
    !case2: h=1 and nonbinding cap         
    try = 2
    IF (choice_fe.EQ.1) try = 3

    hval  = 1d0
    p0    = wbar/alpha
    ctval = ct_global
    cnval = (((1d0-omegac)/omegac/p0)**(1d0/(1d0+mu)))*ctval

    IF (taxrate.GT.1d0) THEN
        IF (chi.EQ.0d0) THEN
            try=4
        ELSE
            try=5
        ENDIF
    ENDIF

    GO TO 1013
    
ELSEIF ((try.GT.1).AND.(try.LT.4)) THEN  
      
    IF (objective_u_unemp.GT.objective_u_unempold) THEN
        objective_u_unempold= objective_u_unemp
        
        if (choice_glob.EQ.1) THEN
            gn_global2old = gn_global2    
            x_global2old  = x_global2 
            ct_global2old = ct_global2
        ENDIF
    ENDIF

    IF (try.EQ.2) THEN

        !case 3: h<1 and binding tax cap
        !new with lumpsum taxes
        IF (taxrate.GT.1d0) THEN
            try =100
        ELSE
            IF (i_num_star2.EQ.num) try = 3
        ENDIF
        
        maxc = ((1d0-omegac)/omegac*alpha/wbar*(ctval**(1d0+mu)))**(alpha/(1d0+alpha*mu))       !b7        
        
        !h < 1
        minc = (alpha/wbar*(1d0-omegac)/omegac)**(1d0/(1d0+mu))*ctval                           !b8
 !       IF (choice_nowbar.EQ.1) THEN
 !           minc = 1d-3
 !           maxc = 1d0
 !        ENDIF
        
        ftype = 302
        IF (i_num_star2.EQ.0)  old_value = AFunction(ftype,a0,aux0,ct_global,minc)
        old_valuemax = 10d+33 

        i_num_star2=i_num_star2
        index_opt = 0
        index_optmax = 0
        difc = (maxc-minc)/(REAL(num)-1d0)

        343 DO i_num=1+i_num_star2,num
            ftype = 302
            ctemp(i_num) = minc+(REAL(i_num)-1d0)*difc     
            new_value = AFunction(ftype,a0,aux0,ct_global,ctemp(i_num))

            sign_new = SIGN(1d0,new_value)
            sign_old = SIGN(1d0,old_value)
            sign_two = SIGN(1d0,sign_new*sign_old)

            if (DABS(new_value) <= DABS(old_valuemax)) then
                index_optmax = i_num
                old_valuemax = new_value        
            end if      
        
            i_num_star2 = i_num
            if (sign_two.EQ.-1d0) then
                index_opt = i_num
                old_value = new_value
                !EXIT
            GOTO 444
            ENDIF
            old_value = new_value
        ENDDO       
 
        444 IF (index_opt.EQ.0) THEN
            bound1=ctemp(MAX(index_optmax-1,1))
            bound2=ctemp(MIN(index_optmax+1,num))
        ELSE
            bound1=ctemp(MAX(index_opt-1,1))
            bound2=ctemp(MIN(index_opt,num))        
        ENDIF

        cnval=BrentRoots(ftype,a0,aux0,ct_global,bound1,bound2,1d-10,1000,fval,niter2,error)  

        IF ((error.EQ.1)) THEN 
            cnval = -10d0
            gnval = -10d0

            IF (i_num_star2.LT.num)  go to 343
            GOTO 1015
        ELSEIF (error.EQ.2) THEN
            PRINT*,'Problem: More Iterations!'
            PAUSE
            STOP
        ENDIF
        
        p0    = (1d0-omegac)/omegac*(ct_global/cnval)**(1d0+mu)
        hval  = (alpha/wbar*p0)**(1d0/(1d0-alpha))

        gnval = hval**alpha-cnval
        tauval= taxrate*(a0+in*p0*hval**alpha)
        
        GO TO 1015     
    ELSEIF (try.EQ.3) THEN         
        !case 4: h=1 and binding tax cap
        IF (chi.EQ.0d0) THEN
            try = 4
        ELSE
            try=5
        ENDIF
                
        hval  = 1d0
        p0    = wbar/alpha
        temp1 = (p0*omegac/(1d0-omegac))**(1d0/(1d0+mu))

        cnval = (a0-taxrate*(a0+in*p0)+p0)/(temp1+p0)
        ctval = temp1*cnval
        gnval = 1d0-cnval
        tauval=p0*gnval-aux0
        taxbar= taxrate*(a0+in*p0)
        
        netbc = tauval-taxbar 

        IF (DABS(netbc).GT.1d-6) tauval=tauval+10d0     !need to impose govt bc with equality
         
        GO TO 1015
    ENDIF
ELSEIF (try.EQ.4) THEN
    IF (objective_u_unemp.GT.objective_u_unempold) THEN
        objective_u_unempold= objective_u_unemp
        
        if (choice_glob.EQ.1) THEN
            gn_global2old = gn_global2    
            x_global2old  = x_global2 
            ct_global2old = ct_global2
        ENDIF
    ENDIF
    ctval=ct_global
    hval = (alpha/wbar*(1d0-omegac)/omegac*(ctval**(1d0+mu)))**(1d0/(1d0+mu*alpha))

    IF ((hval.LE.0d0).OR.(hval.GT.1d0)) THEN
        hval=0d0
        cnval=0d0
        p0=0d0
    ELSE
        cnval=hval**alpha
        p0 = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    ENDIF

    try=5
    gnval  = 0d0
    tauval = -aux0 

    GO TO 1015
ELSE
    IF (objective_u_unempold.GT.objective_u_unemp) THEN
        objective_u_unemp= objective_u_unempold

        if (choice_glob.EQ.1) THEN
            gn_global2 = gn_global2old    
            x_global2  = x_global2old 
            ct_global2 = ct_global2old 
        ENDIF
    ENDIF
ENDIF

objective_u_unemp_SPAIN = objective_u_unemp
!close(10)


end function objective_u_unemp_SPAIN
!******************************************************************************************
!Compute the objective function in the Belman equation
DOUBLE PRECISION function objective_fun_SPAIN(b,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::b,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
INTEGER :: other_type, i,j,h, t, d, NINTV
REAL(8)::u_fun,q_fun,acum,exp_v_next,value_next,a_next,a_fun,scalar,DCSVAL,bb
REAL(8)::a_threshold,a_max1,a_min1,a_next1,acum2,acum1,g,q,DNORDF,v0_fun,vaut_fun,value_next_tirar
REAL(8)::acum_int,acum_int_lo,acum_int_hi,borrowing, interpolate,cdf,cdf1,cdf_threshold, DNORIN,objective_u_unemp_SPAIN
REAL(8)::exp_v_next_aut,acum3,objective_u_fe_SPAIN
EXTERNAL u_fun,q_fun,DNORDF,a_fun,interpolate,v0_fun,vaut_fun,DCSVAL,DNORIN,udef,objective_u_unemp_SPAIN,objective_u_fe_SPAIN

REAL(8)::value1,value2 ,u_value,ucost,udef

IF (psi1.LT.500d0) THEN
    a_threshold   = a_fun(b,a_initial2)
ELSE
    a_threshold   = -1d4
ENDIF

cdf_threshold = DNORDF((a_threshold - rhoa*a_initial2 - (1-rhoa)* mua)/sigmaa)
bb = b*(1d0-defaultgrid(i_default_global2))

exp_v_next = 0d0
exp_v_next_aut = 0d0
scalar = 1d0 / SQRT(pi_number)

a_min1 = (1-rhoa)*mua + rhoa* a_initial2 - 4*sigmaa
a_max1 = (1-rhoa)*mua + rhoa* a_initial2 + 4*sigmaa

NINTV = ncdf - 1
acum1=0d+0
acum2=0d+0

if (a_threshold < a_min1) then

    do i=1,nquad_v
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = v0_fun(bb,a_next)
        acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
    end do

elseif(a_threshold > a_max1) then

    do i=1,neps
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = vaut_fun(a_next)
        acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
    end do

else
    do i=1,neps

    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
    a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
    value_next = vaut_fun(a_next)
    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next

!        compute integral for a > a_threshold
    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
    a_next1 = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf1) * sigmaa !DCSVAL(cdf1, NINTV, BREAK_eps, CSCOEF_eps)
    value_next = v0_fun(bb,a_next1)
    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next

    end do
end if

exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

IF (i_default_global2.EQ.2) THEN
    acum3=0d+0
        do i=1,neps
            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
            a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
            value_next = vaut_fun(a_next)
            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
        end do
    
        exp_v_next_aut = acum3/(cdfmax - cdfmin)
        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
ENDIF

ucost = defaultgrid(i_default_global2)*udef(a_initial2)

value1 = objective_u_fe_SPAIN(bb,a_initial2,b_initial2,i_default_global2)
value2 = objective_u_unemp_SPAIN(bb,a_initial2,b_initial2,i_default_global2)

u_value = MAX(value1,value2)

objective_fun_SPAIN = -(u_value - ucost) - beta * exp_v_next

if (choice_glob.EQ.1) THEN
    IF (value1.GE.value2) THEN
        reg_global = 1
        x_global    = x_global1
        gn_global    = gn_global1
        ctrad_global = ct_global1  
    ELSE
        reg_global = 2
        x_global     = x_global2
        gn_global    = gn_global2
        ctrad_global = ct_global2  
    ENDIF    
ENDIF

end FUNCTION objective_fun_SPAIN
!********************************************************************************************
subroutine optimize_SPAIN(b_next,v_value,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

implicit none

REAL(8),INTENT(IN)::a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
integer :: MAXFN, i,t, d, index_opt, index_vector(1), num, num_tirar
parameter (num=100, num_tirar = 500)        !300
DOUBLE PRECISION :: b_next, v_value, objective_fun_SPAIN, new_value, q_fun, old_value, b, u_fun,q,vector(nb),&
                    b_next_grid(num), STEP, BOUND, XACC, b_next_guess
real(8)::b_next_inf,b_next_sup 
EXTERNAL objective_fun_SPAIN, q_fun, u_fun, DUVMIF

REAL(8)::objective_u_unemp,objective_u_fe,value1,value2
EXTERNAL::objective_u_unemp,objective_u_fe

IF (psi1.GT.500d0) THEN
    b_next_inf = MAX(bmin,b_initial2-0.2d0) !bmin !b_initial2-0.25d0
    b_next_sup = MIN(bmax,b_initial2+0.2d0)
ELSE
    b_next_inf = bmin 
    b_next_sup = bmax
ENDIF

old_value = 10d+33
do i=1,num

   b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
    new_value = objective_fun2(b_next_grid(i))
   
   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
  end if
end do

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value
ELSE
    call DUVMIF(objective_fun2, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
    v_value = -objective_fun2(b_next)
    !value1 = objective_u_fe(b_next,a_initial2,b_initial2,i_default_global2)
    !value2 = objective_u_unemp(b_next,a_initial2,b_initial2,i_default_global2)
end if

CONTAINS

REAL(8) FUNCTION objective_fun2(x)

REAL(8),INTENT(IN)::x

objective_fun2 = objective_fun_SPAIN(x,a_initial2,b_initial2,i_default_global2)

END FUNCTION
100 end subroutine
!************************************************************************************************************
SUBROUTINE pol_fcns_SPAIN
USE GLOBAL

IMPLICIT NONE

REAL(8)::a_initial,b_initial,a0,aux0
INTEGER::i_default_global
REAL(8)::a,b
INTEGER::i_a,i_b,i_bp,i,j,i0(1)
REAL(8)::wmax,p0max,h0max,gn0max,ct0max,tau0max,ww0max,q_fun,objective_fun,interpolate,v_valor,b_next
REAL(8)::ctval,cnval,p0,q,ca1(na2,nb2),q1_nodef(na2,nb2)
EXTERNAL::q_fun,objective_fun,interpolate
REAL(8)::at0
INTEGER::i_at0_spain(9),i_spain,year_spain
CHARACTER(100)::fileplace


do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0
choice_glob=1

OPEN(1,FILE='output//bgrid2.txt')           
DO i = 1,nb2
    WRITE(1,*) bgrid2(i)
ENDDO
CLOSE(1)

OPEN(1,FILE='output//agrid2.txt')           
DO i = 1,na2
    WRITE(1,*) agrid2(i)
ENDDO    
CLOSE(1)

CALL BUILD_TRANSITION(na2,agrid2,rhoa,sigmaa,mua,prob)
   
OPEN(99,FILE='output//gdp_data_spain2007-2013.txt')     !option (B): through 2013

DO i=1,9
    READ(99,*) at0
    i_at0_spain(i) = LOCATE2(agrid2,at0)
ENDDO
CLOSE(99)

endowment_spain(1:5)=0d0
endowment_spain(6)=0.015d0*4.917d0
endowment_spain(7)=0.015d0*4.917d0
endowment_spain(8:9)=0.0d0

DO i_spain=6,9
    year_spain = 2006+i_spain

    WRITE(fileplace,'("output//",I4)') year_spain

    DO i_a=i_at0_spain(i_spain),i_at0_spain(i_spain)
    
    a_initial = agrid2(i_a)
    a0 = DEXP(a_initial)

    endowment = 0d0
    b_initial = 0d0 
    i_default_global = 2  !default occurs
    b = 0d0

    v_valor = -objective_fun(b,a_initial,b_initial,i_default_global)

    vaut1(i_a)  = v_valor
    ctval  = ctrad_global 
    cnval  = x_global**alpha-gn_global
    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    
    paut1(i_a)  = p0
    haut1(i_a)  = x_global
    gnaut1(i_a) = gn_global
    ctaut1(i_a) = ctval
    tauaut1(i_a)= taxrate*(a0+in*p0*x_global**alpha)
    waut1(i_a)  = alpha*p0*x_global**(alpha-1d0)
    
    DO i_b=1,nb2    
            
    i_default_global = 1
    b_initial = bgrid2(i_b)
    endowment = endowment_spain(i_spain)
    
    b = interpolate(b_initial, a_initial, b_next_matrix)
    i_bp = LOCATE2(bgrid2,b)
    
    IF (welfare_case.EQ.1) THEN
        b_next = (1d0-lambda)*b_initial
        v_valor= -objective_fun(b_next,a_initial,b_initial,i_default_global)
    ELSE
        call optimize_SPAIN(b_next, v_valor,a_initial,b_initial,i_default_global)
    ENDIF
    
    vcon1(i_a,i_b)  = v_valor
    bnext1(i_a,i_b) = b_next

    q = q_fun(b_next,a_initial)
    aux0 = -q*(b_next-(1d0-lambda)*b_initial ) + payoff*b_initial            !p*gn - T

    ctval  = ctrad_global
    cnval  = x_global**alpha-gn_global
    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)    
    bp1(i_a,i_b)    = LOCATE2(bgrid2,b_next)
        
    q1_nodef(i_a,i_b) = q_fun(b_next,a_initial)
    p1(i_a,i_b)     = p0
    h1(i_a,i_b)     = x_global
    gn1(i_a,i_b)    = gn_global
    ct1(i_a,i_b)    = ctval
    q1(i_a,i_b)     = q_fun(b_initial,a_initial)
    ca1(i_a,i_b)    = aux0
    tau1(i_a,i_b)   = taxrate*(a0+in*p0*x_global**alpha)
    w1(i_a,i_b)     = alpha*p0*x_global**(alpha-1d0)                    

    def1(i_a,i_b) = 1d0
    IF (vaut1(i_a).GT.vcon1(i_a,i_b)) def1(i_a,i_b) = 0d0
ENDDO
ENDDO


OPEN(1, FILE=TRIM(ADJUSTL(fileplace))//'//gvcon.txt')
OPEN(2, FILE=TRIM(ADJUSTL(fileplace))//'//gp.txt')
OPEN(3, FILE=TRIM(ADJUSTL(fileplace))//'//gh.txt')
OPEN(4, FILE=TRIM(ADJUSTL(fileplace))//'//gibp.txt')
OPEN(5, FILE=TRIM(ADJUSTL(fileplace))//'//gct.txt')
OPEN(6, FILE=TRIM(ADJUSTL(fileplace))//'//gbnext.txt')
OPEN(7, FILE=TRIM(ADJUSTL(fileplace))//'//ggn.txt')
OPEN(8, FILE=TRIM(ADJUSTL(fileplace))//'//gtau.txt')    
OPEN(9, FILE=TRIM(ADJUSTL(fileplace))//'//gw.txt')    
OPEN(10,FILE=TRIM(ADJUSTL(fileplace))//'//gq.txt')    
OPEN(11,FILE=TRIM(ADJUSTL(fileplace))//'//gdef1.txt')    
OPEN(12,FILE=TRIM(ADJUSTL(fileplace))//'//gca.txt')    
OPEN(13,FILE=TRIM(ADJUSTL(fileplace))//'//gq_nodef.txt')    
OPEN(14,FILE=TRIM(ADJUSTL(fileplace))//'//bgrid2.txt')    

DO j = 1,nb2
    WRITE(14,*) bgrid2(j)
ENDDO

CLOSE(14)
DO i = i_at0_spain(i_spain),i_at0_spain(i_spain) !1,na2
DO j = 1,nb2
    WRITE(1,*) vcon1(i,j)
    WRITE(2,*) p1(i,j)
    WRITE(3,*) h1(i,j)
    WRITE(4,'(I7)') bp1(i,j)
    WRITE(6,*) bnext1(i,j)
    WRITE(5,*) ct1(i,j)
    WRITE(7,*) gn1(i,j)
    WRITE(8,*) tau1(i,j)
    WRITE(9,*) w1(i,j)
    WRITE(10,*) q1(i,j)
    WRITE(11,*) def1(i,j)
    WRITE(12,*) ca1(i,j)
    WRITE(13,*) q1_nodef(i,j)
ENDDO
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(4)
CLOSE(5)
CLOSE(6)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)

OPEN(1,FILE=TRIM(ADJUSTL(fileplace))//'//gvaut.txt')
OPEN(2, FILE=TRIM(ADJUSTL(fileplace))//'//gpaut.txt')
OPEN(3, FILE=TRIM(ADJUSTL(fileplace))//'//ghaut.txt')
OPEN(5, FILE=TRIM(ADJUSTL(fileplace))//'//gctaut.txt')
OPEN(7, FILE=TRIM(ADJUSTL(fileplace))//'//ggnaut.txt')
OPEN(8, FILE=TRIM(ADJUSTL(fileplace))//'//gtauaut.txt')    
OPEN(9, FILE=TRIM(ADJUSTL(fileplace))//'//gwaut.txt')    
            
DO i = i_at0_spain(i_spain),i_at0_spain(i_spain) !1,na2
    WRITE(1,*) vaut1(i)
    WRITE(2,*) paut1(i)
    WRITE(3,*) haut1(i)
    WRITE(5,*) ctaut1(i)
    WRITE(7,*) gnaut1(i)
    WRITE(8,*) tauaut1(i)
    WRITE(9,*) waut1(i)
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(5)
CLOSE(7)
CLOSE(8)
CLOSE(9)

PRINT*,'COMPLETED POL FCNS FOR YEAR',2006+i_spain
    
ENDDO

END SUBROUTINE pol_fcns_SPAIN
!****************************************************************************************
